Analisis multivariante 1
Cristopher Valenzuela
27 de febrero de 2022
En este ejercicio vamos a estudiar la tasa de enfermedades transmitidas por los alimentos causadas por el consumo de hortalizas y frutas en la Unión Europea y Estados Unidos con los datos del trabajo de Callejón et al.[1] (2015).
Las enfermedades del análisis fueron Norovirus, Salmonella spp, Escherichia coli, Campylobacter spp, Shigella spp, Clostridium spp, Staphylococcus spp, Yersinia spp, Bacillus spp, otros virus y otros microorganismos. Las verduras eran Ensalada (todos los productos relacionados con la ensalada), Hojas (todos los productos relacionados con las hojas), Tomate y Otras verduras. Las frutas fueron Germinados (todos los productos relacionados con los germinados), Bayas, Melón, Zumos y Otras frutas. Los datos se encuentran en el archivo vegetables.zip. Aunque se pueden cargar desde el archivo vegetables.mat de Matlab, se recomienda leerlos desde los otros archivos. El archivo vegetables.dat contiene tres columnas de datos. La primera columna es la variable illness (enfermedad), la segunda es la variable veg.fruit (verduras/frutas). Y, por último, la tercera columna es la variable region (región). Los nombres de todos los elementos que necesitamos se hallan en el archivo vegetables_labels.txt de forma ordenada según la codificación.
- Obtener las tablas de contingencia entre enfermedades y verduras/frutas separadas por región. brevemente.
#cargamos los datos
vegetables <- read.table("vegetables.dat")
noms <- read.table("vegetables_labels.txt")
#eliminamos los datos en blanco
rm.blanks <- function(x) gsub(" ","",x)
noms <- apply(noms,1,rm.blanks)
colnames(vegetables) <- c("illness","veg.fruit","region")
vegetables$illness <- factor(vegetables$illnes, labels=noms[1:11])
vegetables$veg.fruit <- factor(vegetables$veg.fruit, labels=noms[12:20])
vegetables$region <- factor(vegetables$region, labels=noms[21:22])
#hacemos tablas de contingencia para cada region.
attach(vegetables)
tabla.EU <- table(illness[region=="EU"],veg.fruit[region=="EU"])
tabla.EU##
## Salad Leafy Tomato Other-vegetables Sprouts Berries
## Norovirus 15 26 1 9 0 55
## Salmonella 8 12 1 3 11 0
## E-coli 3 0 0 1 3 0
## Campylobacter 2 1 0 0 0 0
## Shigella 1 0 0 3 0 0
## Clostridium 0 1 0 6 0 0
## Staphylococcus 0 0 0 3 1 0
## Yersinia 0 0 0 3 0 0
## Bacillus 2 0 0 3 0 1
## Other-viruses 2 0 0 0 2 0
## Other-microorganisms 0 0 0 9 0 0
##
## Melon Juices Other-fruits
## Norovirus 0 0 2
## Salmonella 1 1 3
## E-coli 0 0 0
## Campylobacter 0 0 0
## Shigella 0 0 1
## Clostridium 0 0 0
## Staphylococcus 0 0 0
## Yersinia 0 0 0
## Bacillus 0 0 0
## Other-viruses 0 0 1
## Other-microorganisms 0 0 0Tabla de contingencia para USA
tabla.USA <- table(illness[region=="USA"],veg.fruit[region=="USA"])
tabla.USA##
## Salad Leafy Tomato Other-vegetables Sprouts Berries
## Norovirus 97 62 5 9 0 5
## Salmonella 8 8 17 3 14 2
## E-coli 10 22 0 0 4 2
## Campylobacter 4 2 1 0 0 0
## Shigella 1 2 0 0 0 0
## Clostridium 0 0 0 0 0 0
## Staphylococcus 2 0 0 0 0 0
## Yersinia 0 0 0 0 0 0
## Bacillus 1 0 0 0 0 0
## Other-viruses 0 1 1 1 0 0
## Other-microorganisms 0 0 0 0 2 0
##
## Melon Juices Other-fruits
## Norovirus 9 3 33
## Salmonella 14 0 5
## E-coli 0 6 2
## Campylobacter 1 0 1
## Shigella 0 0 0
## Clostridium 0 0 0
## Staphylococcus 0 0 0
## Yersinia 0 0 0
## Bacillus 0 0 1
## Other-viruses 0 0 2
## Other-microorganisms 1 0 0- Realizar un análisis de correspondencias1 de las tablas, representar los resultados y comentarlos brevemente.
library(ca)
ca.EU <- ca(tabla.EU)
round(ca.EU$sv^2, 3) # valores propios## [1] 0.501 0.425 0.069 0.042 0.035 0.003 0.000 0.000plot(ca.EU)
Tambien realizamos un plot asimetrico
plot(ca.EU, map = "rowprincipal")
Del mismo modo se puede hacer el otro gráfico asimétrico:
plot(ca.EU, map = "colprincipal")
Realizamos los mismo graficos para USA, pero antes eliminamos datos faltantes
ca.USA <- ca(tabla.USA[-c(6,8),])
round(ca.USA$sv^2, 3) # valores propios## [1] 0.384 0.139 0.041 0.022 0.008 0.003 0.001 0.000plot(ca.USA)
Tambien hacemos graficos asimetricos
plot(ca.USA, map = "rowprincipal")
plot(ca.USA, map = "colprincipal")
Para USA los dos primeros componentes explican el 87.2% de los datos, mientras que para la UE explican el 86.1%.
El Norovirus se muestra como responsable de la mayoría de los brotes relacionados con productos agrícolas, seguido de la Salmonela. El Norovirus está relacionado principalmente con el consumo de ensalada en los Estados Unidos y de bayas en la Unión Europea. La Salmonella fue la principal causa de brotes de productos agrícolas en varios estados de Estados Unidos y fue el patógeno implicado en la mayoría de los brotes asociados al consumo de coles, tomates y melones. Como se refleja en los gráficos, el patrón de los brotes de productos frescos difería en Estados Unidos y en la Unión Europea según el tipo de microorganismo y el vehículo alimentario implicado.
Corradino[2] utilizó el MDS para estudiar la estructura de proximidades en una colonia de monos japoneses. Las observaciones se realizaron en un grupo social de 14 monos japoneses durante un período de un año. Los catorce monos se nombran y describen en la Tabla 1. Se observaron las relaciones de proximidad cada 60 segundos. Si dos monos se encontraban a menos de 1.5 metros de distancia y se toleraban, se decía que estaban “cerca”. Se calcularon las disimilaridades para cada par de monos en función del tiempo que la pareja estuvo cerca la una de la otra. Las disimilaridades se estudiaron por separado en la temporada de cría (monk_85.dis) y fuera de ésta (monk_84.dis).
- Utilizar las disimilaridades del archivo monk_84.dis2 para definir un objeto del tipo dist en R o una matriz simétrica con esas disimilaridades. Observaremos el orden de la tabla 1 que se puede obtener del archivo monkeys.dat. El vector con las 91 disimilaridades debe rellenar la matriz triangular inferior del objeto por columnas. También habrá que poner nombres a las filas y columnas de la matriz.
Lo mismo para el archivo monk_85.dis con las disimilaridades en la temporada de cría.
El primer paso es crear la matriz de distancia para los datos
monkeys <- read.csv("monkeys.dat", sep=";")
monkeys$age <- factor(monkeys$age)
monkeys$sex <- factor(monkeys$sex)
df <- read.table("monk_84.dis", header = FALSE, sep="", skip=2)
disim <- as.numeric(df$V1[1:91])
mat <- matrix(0, ncol=14, nrow=14)
mat[lower.tri(mat)] <- disim
mat <- mat + t(mat)
colnames(mat) <- rownames(mat) <- monkeys$alias
d.84 <- as.dist(mat)
d.84## ALFA FRAN FELL PANC ISA GILD BETI OLGA ORSE ROSS DIVO CIST ELET
## FRAN 126
## FELL 55 211
## PANC 367 60 26
## ISA 55 39 21 12
## GILD 63 39 17 23 1
## BETI 37 74 73 32 11 3
## OLGA 27 33 24 2 28 3 3
## ORSE 96 90 68 8 14 26 24 6
## ROSS 62 30 8 51 20 1 35 8 42
## DIVO 27 18 51 100 10 17 8 5 28 0
## CIST 4 10 100 38 1 3 3 17 39 3 16
## ELET 3 19 38 16 0 35 7 16 4 4 62 6
## EVA 22 5 16 46 8 9 1 6 46 16 9 26 436Realizamos la matriz de distancias para la otra temporada.
df <- read.table("monk_85.dis", header = FALSE, sep="", skip=2)
disim <- as.numeric(df$V1[1:91])
mat <- matrix(0, ncol=14, nrow=14)
mat[lower.tri(mat)] <- disim
mat <- mat + t(mat)
colnames(mat) <- rownames(mat) <- monkeys$alias
d.85 <- as.dist(mat)
d.85## ALFA FRAN FELL PANC ISA GILD BETI OLGA ORSE ROSS DIVO CIST ELET
## FRAN 80
## FELL 44 161
## PANC 156 120 27
## ISA 49 18 10 21
## GILD 109 117 22 4 2
## BETI 77 113 23 40 78 11
## OLGA 83 105 9 9 7 74 18
## ORSE 53 153 105 41 25 47 69 4
## ROSS 33 44 7 11 4 2 14 16 9
## DIVO 20 29 28 75 20 44 60 14 39 0
## CIST 22 22 20 9 4 56 107 28 62 15 11
## ELET 286 7 27 7 2 41 7 14 3 15 109 34
## EVA 59 17 256 18 12 38 0 146 96 20 12 22 204- Comprobar que las disimilaridades en las dos temporadas forman conjuntos no euclídeos. ¿Qué significa esto?
library(ade4)
is.euclid(d.84,print=F)## Warning in is.euclid(d.84, print = F): Zero distance(s)## [1] FALSE#como nos faltan datos chequeamos donde estan
which(d.84==0)## [1] 54 82Reemplazamos estos datos por un valor pequeño cercano al "0"
d.84[which(d.84==0)] <- 0.0001
which(d.85==0)## [1] 70 82d.85[which(d.85==0)] <- 0.0001Comprobamos si la distancia es euclidea
is.euclid(d.84,print=F)## [1] FALSEEn este caso la distancia no es euclidea.
is.euclid(d.85,print=F)## [1] FALSEEn este caso tampoco es euclidea.
Cuando la disimilaridad es euclídea o una medida directa relativa a ella (tal como el coseno), la mejor elección para el MDS es el MDS métrico. Sin embargo, cuando la disimilaridad no es euclídea o incluso no métrica, debemos decidir entre la solución MDS métrica y la no métrica. El MDS no métrico frecuentemente demuestra ser más realista y mucho mejor en la práctica.
- Realizar el MDS más apropiado (clásico, métrico o no métrico) con las disimilaridades de las dos temporadas por separado. ¿Cual es el stress en cada caso?
set.seed(123)
library(MASS)
iso.mds.84 <- isoMDS(d.84)## initial value 62.369024
## final value 62.351900
## convergediso.mds.84$stress## [1] 62.3519iso.mds.85 <- isoMDS(d.85)## initial value 47.386485
## final value 47.345903
## convergediso.mds.85$stress## [1] 47.3459library(vegan)## Loading required package: permute## Loading required package: lattice## This is vegan 2.5-7nmds.84 <- metaMDS(d.84, engine="monoMDS", autotransform=FALSE, trace = FALSE)
nmds.84$stress*100## [1] 20.91301nmds.85 <- metaMDS(d.85, engine="monoMDS", autotransform=FALSE, trace = FALSE)
nmds.85$stress*100## [1] 22.84126Aunque con valores elevados de stress, es evidente que la mejor solución se obtiene con la función metaMDS().
- Dibujar la representación de puntos con las mejores soluciones del apartado anterior, una para cada temporada.
plot(nmds.84$points[,1], nmds.84$points[,2], type = "n",
xlab="NMDS1", ylab="NMDS2")
title("Temporada sin cría")
text(nmds.84$points[,1], nmds.84$points[,2],
col=ifelse(monkeys$sex=="male","blue","red"),
labels = rownames(nmds.84$points), cex=0.5)
plot(nmds.84$points[,1], nmds.84$points[,2], type = "n",
xlab="NMDS1", ylab="NMDS2")
title("Temporada de cría")
text(nmds.85$points[,1], nmds.85$points[,2],
col=ifelse(monkeys$sex=="male","blue","red"),
labels = rownames(nmds.85$points), cex=0.5)
Observamos una mayor proximidad en la representación de la temporada de cría. Los machos aparecen agrupados en las dos configuraciones, aunque con alguna excepción. Las hembras también forman un grupo, es decir, se relacionan más entre ellas en las dos temporadas.
En el trabajo de Corradino se deducía que los machos estaban más cerca en la temporada de cría que fuera de ella. En nuestros gráficos este efecto no aparece.
- Comparar las dos configuraciones de las temporadas con la función procrustes().
¿Cual es el mono con menor residuo entre las dos configuraciones?
proc <- procrustes(nmds.84,nmds.85)
plot(proc)
El menor residuo de la representación es para Gilda.
sort(residuals(proc))[1:3]## GILD EVA DIVO
## 31.67155 41.72386 61.20314- Como el stress es alto, podemos optar por una representación en tres dimensiones.
En primer lugar procedemos a calcular el MDS no métrico con 3 coordenadas.
nmds.85.3 <- metaMDS(d.85, engine="monoMDS", k = 3, autotransform=FALSE, trace = FALSE)
nmds.85.3$stress*100## [1] 14.81673El stress es más bajo que con dos dimensiones.
library(plotly)## Loading required package: ggplot2##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:MASS':
##
## select## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutdatos <- data.frame(nmds.85.3$points)
fig3D <- plot_ly(data=datos, x = ~MDS1, y = ~MDS2, z = ~MDS3,
color = ~monkeys$sex, colors = c("red","blue"))
fig3D## No trace type specified:
## Based on info supplied, a 'scatter3d' trace seems appropriate.
## Read more about this trace type -> https://plotly.com/r/reference/#scatter3d## No scatter3d mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plotly.com/r/reference/#scatter-modeCon la rotación adecuada se puede ver la separación total entre machos y hembras.
Normalmente, el análisis de componentes principales se lleva a cabo calculando los valores propios y los vectores propios de la matriz de correlaciones. Con n casos y p variables, si escribimos X para la matriz n×p estandarizada para que las columnas tengan media cero y desviación estándar unitaria, encontramos los valores propios y los vectores propios de la matriz p × p X0X (que es n o n − 1 veces la matriz de correlaciones dependiendo de qué denominador se haya utilizado al calcular las desviaciones estándar). Con los datos de los gorriones del ejercicio 3 de los ejercicios propuestos para PCA, comprobar las siguientes afirmaciones:
- El primer vector propio da las cargas (loadings) de cada variable en la primera componente principal, el segundo vector propio da las cargas en la segunda componente y así sucesivamente. Comprobar este resultado con la función princomp().
load("gorriones.RData")
X <- gorriones[,1:5]
X <- scale(X)
eig <- eigen(crossprod(X), symmetric = TRUE)
pca <- princomp(X)
pca$loadings##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## x1 0.452 0.690 0.420 0.374
## x2 0.462 -0.300 0.341 -0.548 -0.530
## x3 0.451 -0.325 -0.454 0.606 -0.343
## x4 0.471 -0.185 -0.411 -0.388 0.652
## x5 0.398 0.876 -0.178 -0.192
##
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## SS loadings 1.0 1.0 1.0 1.0 1.0
## Proportion Var 0.2 0.2 0.2 0.2 0.2
## Cumulative Var 0.2 0.4 0.6 0.8 1.0 eig$vectors## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.4517989 0.05072137 0.6904702 0.42041399 -0.3739091
## [2,] 0.4616809 -0.29956355 0.3405484 -0.54786307 0.5300805
## [3,] 0.4505416 -0.32457242 -0.4544927 0.60629605 0.3427923
## [4,] 0.4707389 -0.18468403 -0.4109350 -0.38827811 -0.6516665
## [5,] 0.3976754 0.87648935 -0.1784558 -0.06887199 0.1924341- Escribiendo las cargas de las primeras q componentes como columnas de la matriz p×q B, la matriz n × q P de las puntuaciones (scores) de las componentes principales de los individuos se obtiene aplicando las cargas a la matriz de datos original, es decir, P = XB.
scores <- X %*% eig$vectors
head(scores)## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.06428901 -0.60083713 -0.1712334 -0.515825561 0.5487904
## [2,] -2.18031283 -0.44230082 0.4000696 -0.645459959 0.2310766
## [3,] -1.14556567 0.01925412 -0.6761269 -0.716298164 0.2088714
## [4,] -2.31106565 0.17199267 -0.3059621 0.149289289 0.4781034
## [5,] -0.29504203 -0.66520783 -0.4742138 -0.545862110 0.2444780
## [6,] 1.91626198 -0.59525444 0.6209330 0.006608669 -0.2855166head(pca$scores)## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## [1,] 0.06428901 -0.60083713 -0.1712334 -0.515825561 -0.5487904
## [2,] -2.18031283 -0.44230082 0.4000696 -0.645459959 -0.2310766
## [3,] -1.14556567 0.01925412 -0.6761269 -0.716298164 -0.2088714
## [4,] -2.31106565 0.17199267 -0.3059621 0.149289289 -0.4781034
## [5,] -0.29504203 -0.66520783 -0.4742138 -0.545862110 -0.2444780
## [6,] 1.91626198 -0.59525444 0.6209330 0.006608669 0.2855166- La matriz de suma de cuadrados y productos, P'P = D, es diagonal con elementos iguales a los primeros q valores propios de de la matriz X0X, de modo que las varianzas de las componentes principales pueden obtenerse dividiendo los valores propios por n o n − 1.
D <- crossprod(scores)
eig$values## [1] 173.566961 25.512196 18.548378 14.475145 7.897321diag(D)## [1] 173.566961 25.512196 18.548378 14.475145 7.897321diag(D)/49## [1] 3.5421829 0.5206571 0.3785383 0.2954111 0.1611698pca$sdev^2## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## 3.5421829 0.5206571 0.3785383 0.2954111 0.1611698- Sustituimos los valores desconocidos (missings) de la matriz X (la de los controles) por la media de cada columna. El siguiente es estandarizar los datos para que tengan media cero y desviación estándar uno y finalmente calcular la matriz XX'.
load("SNP.RData")
controls <- rownames(subject.support)[subject.support$cc==0]
pop <- subject.support[controls,"stratum"]
pop.all <- subject.support[,"stratum"]
use <- seq(1, ncol(snps.10), 10)## Loading required package: snpStats## Loading required package: survival## Loading required package: Matrixctl10 <- snps.10[controls, use]
all10 <- snps.10[,use]
X <- as(ctl10, "numeric")
X.all <- as(all10,"numeric")
dim(X)## [1] 500 2851La primera dificultad es substituir los valores missing por la media.
for(i in 1:ncol(X)){
X[is.na(X[,i]), i] <- mean(X[,i], na.rm = TRUE)
}Estandarizamos los datos
X <- scale(X)Calculamos la matriz XX'
XXt <- tcrossprod(X)
head(XXt[1:10,1:10])## jpt.869 jpt.862 jpt.948 ceu.564 ceu.904 ceu.665 jpt.663
## jpt.869 2601.7764 138.9807 192.9046 -283.5339 -414.6340 -235.4016 247.0385
## jpt.862 138.9807 2775.1951 164.8674 -218.8962 -306.4598 -280.7771 202.4891
## jpt.948 192.9046 164.8674 2531.5733 -361.6866 -311.3615 -278.4234 309.6210
## ceu.564 -283.5339 -218.8962 -361.6866 2908.3334 439.0646 329.4606 -162.9436
## ceu.904 -414.6340 -306.4598 -311.3615 439.0646 2935.7316 313.8808 -296.6975
## ceu.665 -235.4016 -280.7771 -278.4234 329.4606 313.8808 2836.9213 -326.4090
## ceu.977 jpt.637 ceu.897
## jpt.869 -311.0834 329.9604 -95.97765
## jpt.862 -345.0864 304.8939 -306.84784
## jpt.948 -284.2715 316.8407 -261.73112
## ceu.564 260.8681 -310.9871 160.09227
## ceu.904 229.0223 -403.5789 308.65566
## ceu.665 355.8414 -296.9524 329.36431- Calcular los valores y vectores propios de la matriz XX' y, a partir de ellos, la matriz U. Comparar gráficamente las puntuaciones de la primera componente para los controles según la población.
Hacer el mismo gráfico con la segunda y tercera componentes. ¿Qué componente separa mejor las poblaciones?
eig <- eigen(XXt, symmetric = TRUE)
eig$values[1:5]## [1] 146789.918 8589.204 8380.229 8236.379 8111.603U <- eig$vectorsoldpar <- par(mfrow=c(1,3))
boxplot(U[,1] ~ pop, ylab="PC1")
boxplot(U[,2] ~ pop, ylab="PC2")
boxplot(U[,3] ~ pop, ylab="PC3")
par(oldpar)- Calcular las cargas B según se ha explicado antes.
eig$values[eig$values<0] <- 0
v <- 1/sqrt(eig$values)
m <- crossprod(X, U)
# B <- m %*% diag(v) # poco eficiente
# B <- t(t(m) * v) # más eficiente
B <- m * rep(v, rep(nrow(m), ncol(m))) # mejor
B[1:10,1:10]## [,1] [,2] [,3] [,4] [,5]
## rs7909677 0.0003252226 -9.331756e-03 0.015445805 0.003947906 -0.011546724
## rs4880781 -0.0055876444 7.507270e-03 -0.000543355 -0.008346595 -0.003828963
## rs2448365 -0.0076646829 1.733169e-02 -0.008138940 0.007698970 0.007641585
## rs7919436 0.0209014311 1.102595e-02 -0.025856734 0.013353554 0.024974745
## rs12357593 0.0301087215 1.708045e-03 -0.022693905 0.014494197 0.022537829
## rs10904173 -0.0032615624 1.569638e-02 0.017220464 0.002238865 -0.014006605
## rs11252693 0.0107882727 -4.632472e-05 -0.012653201 0.011259890 0.006176045
## rs11253096 0.0259380458 -1.762775e-02 0.006504495 0.010741717 -0.010751181
## rs17159711 0.0150469759 1.284841e-02 0.010802802 0.007726009 0.030637888
## rs4881518 0.0168951870 -1.775938e-02 -0.008888053 0.018592907 0.004927426
## [,6] [,7] [,8] [,9] [,10]
## rs7909677 0.009234650 0.009592387 -0.0006004881 0.009523450 -0.014859529
## rs4880781 -0.025199111 0.001262191 0.0136031063 -0.034072213 -0.011145068
## rs2448365 -0.009577234 0.025332822 0.0176713475 0.005336778 -0.009796594
## rs7919436 0.004708050 -0.017783670 0.0129060916 0.035684897 -0.009938182
## rs12357593 0.017645674 -0.035232563 0.0028693908 0.036582396 -0.002629847
## rs10904173 0.008486269 0.013840947 0.0097519493 -0.023745230 -0.017696460
## rs11252693 0.019020146 -0.019180868 -0.0004079738 0.033007216 0.014243545
## rs11253096 0.004788495 -0.014961988 0.0022505070 -0.024822305 -0.008882354
## rs17159711 -0.011462165 0.022790028 -0.0054857879 0.008332284 -0.009762654
## rs4881518 0.022255034 -0.006821502 0.0213503830 0.025230784 0.001849588- Calcular las puntuaciones de las dos primeras componentes de todas las observaciones con las cargas de los controles.
Se controlan los valores faltantes y se estandarizan los datos
for(i in 1:ncol(X.all)){
X.all[is.na(X.all[,i]), i] <- mean(X.all[,i], na.rm = TRUE)
}
X.all <- scale(X.all)Calculamos las puntuaciones y hacemos el grafico
scores.all.12 <- X.all %*% B[,1:2]
plot(scores.all.12, xlab="PC1", ylab="PC2",
col=ifelse(pop.all=="CEU","red","blue"))
Las dos poblaciones quedan también separadas por la primera componente.